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Ill 


Abstract 

The  initialization  problem  is  defined  as  the  problem  of 
obtaining  the  initial  data  that  are  required  in  order  to  solve 
a  well  posed  initial-and-boundary-value  problem  for  the  equa- 
tions of  large  scale  dynamical  meteorology.   In  our  treatment 
of  this  problem  complete  knowledge  at  a  given  instant  of  the 
atmospheric  temperature  and  surface-pressure  fields,  as  well  as 
of  their  time  derivatives,  is  assumed  to  be  at  our  disposal;  and 
the  problem  of  computing  the  velocity  field  at  that  instant  from 
this  knowledge  is  studied. 

The  methods  used  are  based  on  the  exact  mathematical  treat- 
ment of  the  differential  equations  of  atmospheric  motion,  rather 
than  on  perturbation- type  or  numerical  methods.   A  simple  solu- 
tion of  the  initialization  problem  thus  formulated  is  exhibited 
in  a  particular  case.   It  is  shown,  however,  that  in  general  the 
problem  is  indeterminate  and  additional  information  is  required 
to  uniquely  find  the  velocity  field.   Four-dimensional  data 
assimilation  is  also  discussed. 


IV 


1.  Introduction 

Shortly  after  the  development  of  computers  made  the  numeri- 
cal integration  of  the  primitive  equations^  over  a  global  grid 
feasible,  it  became  apparent  that  one  of  the  main  problems  in 
such  an  integration  was  to  provide  the  initial  data  needed  to 
start  it.   Indeed,  the  primitive  equations,  in  x,y, z-coordinates, 
say,  form  a  totally  hyperbolic  system,  when  viscosity  and  heat 
conduction  are  neglected.    A  well-posed  problem  for  this  system 
would  be  to  specify  throughout  the  global  atmosphere  the  velocity 
components  u,  v,  w,  and  the  temperature  T  at  a  given  time  t  =  t^, 
as  well  as  certain  conditions,  into  the  details  of  which  we  will 
not  enter,  at  the  upper  and  at  the  lower  boundary  of  the  atmo- 
sphere.  The  complete  and  sufficiently  accurate  specification  of 
such  initial  conditions  became  known  as  the  initialization 
problem. 

Thus  the  problem  falls  into  two  parts:   (i)  the  complete- 
ness, and  (ii)  the  accuracy  of  the  initial  data.   In  discussing 
the  different  aspects  of  this  problem,  we  will  try  to  distinguish 
between  the  physical  behavior  of  the  atmosphere,  the  mathematical 
(differential)  models  describing  this  behavior,  and  the  numerical 
(difference)  models  used  to  approximate  the  mathematical  model. 
The  behavior  of  a  numerical  model  is  often  related  directly  to 


'  The  Eulerian  equations  of  motion  with  the  vertical  momentum 
equation  replaced  by  the  hydrostatic  ass\Amption. 

Assiming  the  surface  pressure  is  known,  the  pressure  p  can  be 
obtained  directly  by  integrating  the  hydrostatic  equation,  which 
we  exclude  from  this  discussion. 


and  compared  with  that  of  the  atmosphere,  the  Intermediate  level 
of  approximation  provided  by  the  mathematical  model  from  which 
it  was  derived  being  omitted  from  the  discussion  or  the  behavior 
of  the  mathematical  model  being  tacitly  implied  to  be  the  same 
as  either  that  of  the  atmosphere,  or  as  that  of  the  numerical 
model.   We  shall  briefly  mention  first  the  questions  related  to 
the  accuracy  problem. 

When  finite-difference  approximations  to  the  primitive 
equations  started  to  be  used  for  meteorological  research  and 
weather  prediction,  it  was  noted  that  the  numerical  models   simu- 
lated not  only  the  large-scale  motions  of  the  atmosphere,  in 
which  one  was  interested,  but  also  faster  phenomena  of  smaller 
amplitude  which  were  identified  as  the  numerical  counterpart  of 
inertia-gravity  waves,  and  which  are  adequately  represented  for 
instance  by  solutions  of  the  linearized  shallow-fluid  equations 
(Charney,  1955).   These  waves  were  undesirable  and  their  cause 
was  traced  to  the  lack  of  geostrophic  adjustment  in  the  initial 
data  usedo   However,  it  was  soon  realized  that  initial  condi- 
tions, satisfying  any  hitherto  used  balance  equation,  would  not 
prevent  the  existence  and  propagation  of  inertia-gravity  waves 
in  a  primitive-equation  model  (Nitta  and  Hovermale,  1969^  Morel 
et  alo,  1971).   The  mathematical  reason  for  this  is  the  fact 


In  this  context  a  numerical  model  of  the  atmosphere  is 
essentially  a  finite-difference  approximation  to  the  equations 
of  motion.   It  includes  in  general  also  some  terms  which 
represent  the  direct  numerical  parametrization  of  certain 
physical  processes,  but  we  will  not  dwell  upon  this  aspect  of 
the  models,  considering  them  as  difference  schemes  for  the  inte- 
gration of  the  equations  of  motion,  i.e.,  of  the  mathematical 
model. 


that  the  forms  of  the  balance  equation  known  from  previous 
quasi-geostrophic  and  quasi-nondivergent  models  do  not  represent 
an  instantaneous  compatibility  condition  between  wind  and  pres- 
sure field,  consistent  with  the  primitive  equations  describing 
the  motion.   Actually  such  a  balance  equation,  in  its  general 
form  or  in  any  of  its  simplified  forms,  represents  only  an 
approximation  to  such  a  condition,  obtained  by  neglecting  certain 
terms  in  the  equations  from  which  it  is  derived.   Indeed,  the 
atmosphere  itself  is  only  approximately  in  geostrophic  balance, 
rather  than  exactly  in  this  state. 

It  appeared  though  that  the  amplitude  of  the  fast  waves  in 
the  numerical  model  far  exceeded  that  of  the  inertia-gravity 
waves  in  the  atmosphere.   These  large  amplitudes  could  reasonably 
be  attributed  to  an  excessive  imbalance  in  the  initial  conditions, 
due  to  the  use  of  data  incorrectly  measured  or  interpolated. 
Since  this  situation  could  not  be  remedied  by  traditional  synop- 
tic balancing,  Nitta  and  Hovermale  (I969)  proposed  a  dynamical 
balancing.   The  procedure  consisted  essentially  in  running  the 
model  successively  for  short  periods  foorward  and  backward  in  time 
for  a  number  of  cycles.   It  was  expected  that  the  model  would 
simulate  geostrophic  adjustment,  i.e.,  disperse  and  dissipate  the 
fast,  short-scale  waves.   Positive  improvement  was  obtained  with 
this  method. 


For  the  most  general  form  of  the  balance  equation,  as  well  as 
for  some  of  the  more  special  forms  it  takes  when  certain  simpli- 
fying assumptions  are  made,  see  Haltiner  (1971). 


The  reduction  of  the  exaggerated  amplitude  of  small-scale 
motions  produced  by  initial  errors  was  known  to  be  necessary  in 
order  to  prevent  them  from  influencing  too  strongly  in  too  short 
a  time  span  the  large-scale  motions  represented  by  the  n^umerical 
model.   The  relationship  between  the  size  of  initial  errors  and 
the  predictability  range   has  been  first  investigated  systemati- 
cally by  Charney  et  al.  ( I966 ) .   For  the  presently  available 
accuracy  in  measurement,  the  predictability  range  of  the  existing 
numerical  models  of  the  atmosphere  is  believed  to  be  about  two 
weeks,  ass-uming  that  complete  information  with  the  given  accuracy 
were  available. 

This  brings  us  back  to  the  completeness  problem.   In  fact, 
there  is  at  present  no  meteorological  data  gathering  system 
capable  of  supplying  u,  v,  w,  T  over  a  uniform  grid,  like  the 
ones  used  by  general  circulation  models  (GCMs),  at  a  given  time 
(" synoptically" ) .   However,  it  seems  reasonable  to  expect  from 
the  progress  of  satellite  technology  at  least  good  date  coverage 
for  temperatures  in  the  near  future.   Therefore  it  has  been 
suggested  by  Charney  et  al.  ( 1969 )  that  the  wind  field  (u, v,w), 
speaking  vaguely  for  the  moment,  could  be  inferred  from  the 
temperature  field  T  and  some  of  its  past  history. 

The  rationale  for  this  proposal  was  twofold.   First,  the 
experience  with  barotropic  models,  i.e.,  with  systems  of  partial 


The  predictability  range  is  defined  as  the  time  span  after 
which  the  errors  in  the  calculated  state  exceed  a  certain  pre- 
assigned  threshold,  usually  related  to  the  difference  between 
two  "randomly  chosen"  states  of  the  atmosphere. 


differential  equations  in  which  the  tiarotropic  approximation  was 
made,  and  in  which  the  wind  field  could  be  obtained  from  the 
instantaneous  geopotential  field  by  solving  a  time-independent 
balance  equation.   Numerical  integrations  of  these  equations  were 
rather  successful,  approximating  well  certain  aspects  of  the  be- 
havior of  the  atmosphere.   Second,  the  well-known  mathematical 
fact  that,  in  a  linear  hyperbolic  system  of  n  first-order  partial 
differential  equations  with  constant  coefficients,  (n-1)  depend- 
ent variables  can  be  eliminated,  thereby  reducing  the  system  to 
a  single  equation  for  the  n-th  variable.   Then  the  initial  condi- 
tions on  the  (n-1)  eliminated  variables  lead  to  conditions  for 
the  time  derivatives  up  to  order  (n-1)  of  the  n-th  variable 
(e.g.,  Courant  and  Hilbert,  I962,    p.  l4ff.).   In  particular, 
this  is  the  case  in  linear  acoustics,  where  the  velocity  compo- 
nents can  be  eliminated,  to  yield  a  single  equation  for  the  pres- 
sure.  The  discussion  of  these  considerations  will  be  briefly 
taken  up  in  Section  4. 

The  concrete  procedure  proposed  by  Charney  et  al.  (1969) 
was  to  start  the  numerical  integration  of  the  equations  with  the 
"correct"  temperature  data  and  with  velocity  data  approximated  by 
some  other  method,  then,  at  given  time  intervals,  to  replace  the 
calculated  values  of  T  by  the  "correct"  ones  (i.e.,  to  "update" 
T).   In  fact,  however,  no  "correct"  measured  (satellite)  tempera- 
ture data  were  available.   Therefore  the  experiments  were  made  to 
test  the  procedure  by  relying  on  the  concept  of  a  "control  run", 
i.e.,  of  a  numerical  integration  of  the  equations  performed  for  a 
certain  length  of  time,  the  results  of  which  are  referred  to  for 


purposes  of  comparison  as  the  "true  history"  of  the  atmosphere 
over  that  time  span. 

First  the  control  run,  with  certain  initial  values  of  u, 
V,  w,  T,    is  performed  and  the  values  of  T  corresponding  to  this 
run  are  taken  as  "correct".   Then  another  integration  is  started, 
with  the  same  initial  temperature  data,  but  different  initial 
velocity  data,  and  T  is  "updated"  with  values  taken  from  the 
"correct"  temperatures  of  the  control  run.   The  reason  for  up- 
dating the  temperature  is  to  reduce  the  error  in  the  velocity, 
caused  by  using  inaccurate  initial  velocity  data. 

It  was  indeed  observed  that  the  root-mean- square  errors  in 
the  velocity  field  decreased  to  a  certain  non-zero  asymptotic 
value.   However,  in  the  experiments  reported  by  Charney  et  al. 
(1969),  as  well  as  in  similar  experiments  performed  later  (e.g. 
Williamson  and  Kasahara,  1971)  this  value  was  reached  only  after 
a  time  comparable  to  the  predictability  range  of  the  numerical 
models  used  (10  days  to  3  weeks)  and  this  asymptotic  rms  error 
was  non-negligible,  although  considerably  smaller  than  the  rms 
value  of  the  initial  error,  and  a  fortiori  much  smaller  than  the 
rms  error  in  the  wind  field  after  a  similar  period  in  a  run  with- 
out updating.   These  facts  seem  to  limit  the  value  of  the  up- 
dating approach  to  the  completeness  problem  in  spite  of  the 
partial  success  of  the  method. 

More  recently  (Morel  et  al.,  1971^  Mesinger,  1972)  the 
method  of  iterative  dynamic  adjustment,  suggested  by  Nitta  and 
Hovermale  (I969)  in  dealing  with  the  geostrophic  balance  problem, 
was  applied  to  the  data  acquisition  problem.   This  was  done  in 


order  to  circiimvent  the  questions  arising  from  achieving  proper 
initialization  only  at  the  end  of  the  predictability  timeo   Good 
results  were  reported  with  proper  choice  of  finite-difference 
scheme,  period  of  the  back-and-f orth  iteration,  etc.   There  seems 
however  to  be  a  marked  distinction  between  the  results  using  data 
from  a  control  run  of  the  same  model  and  those  when  using  real 
data  or  data  from  a  different  model.   The  results  when  using 
"extraneous"  data  for  iterative  updating  are  much  worse  and 
rather  discouraging  (Morel  et  al.,  1971^  Mesinger,  1972). 

It  may  be  worthwhile  to  stress  that  the  forward-updating 
and  iterative  techniques  discussed  before  were  proposed  when  it 
became  clear  that  traditional  objective  analysis  methods,  which 
had  been  fairly  successful  for  the  "filtered"  (quasi-geostrophic 
and  quadi-nondivergent )  models,  produced  "initialization  shocks" 
and  large-amplitude  short  waves  in  primitive  equation  models. 
However,  the  consensus  seems  to  be  that  up  to  now  these  asynchro- 
nous techniques  are  much  more  computer-time  consuming  and  do  not 
seem  to  have  performed  convincingly  better  than  the  more  sophis- 
ticated ones  among  the  objective-analysis  methods  (including 
least-squares  and  weighted-average  fits  of  observed  data,  varia- 
tional techniques,  weighted  consideration  of  prognostic  computed 
data,  permanence  and  climatology  data,  as  well  as  use  of  geo- 
strophic  and  balanced  data).   Therefore  it  is  still  the  synchro- 
nous objective  analysis  methods  which  are  operational  in  daily 
forecasting  (Haltiner,  1971)  and  it  is  one  of  the  purposes  of 
the  search  for  better  initialization  methods  to  develop  tools 
for  improved  and  extended  operational  forecasts. 


In  view  of  the  difficulties  encountered  by  the  non-synoptic 
techniques,  as  well  as  by  the  traditional  synoptic  techniques, 
we  shall  try  to  pursue  a  different  avenue  of  attack  on  the 
initialization  problem.   As  we  saw,  the  two  aspects  of  the  prob- 
lem, completeness  and  accuracy,  are  strongly  interrelated  and 
any  solution  has  to  take  care  of  both.   This  brings  us  back  to 
the  discussion  of  well-posed  problems  for  the  system  of  equations 
dealt  with,  viz.,  the  primitive  equations  (where  we  assume  the 
surface  pressure  to  be  known  and  hence  ignore  the  hydrostatic 
equation),  i.e.,  to  the  question  of  the  side  conditions  under 
which  the  existence  of  a  unique  solution  of  this  system,  con- 
tinuously depending  on  the  data,  is  guaranteed. 

Given  the  hyperbolic  character  of  the  equations  and  their 
non-linearity  it  seems  mathematically  very  implausible  that  any 
other  problem,  but  the  proper  initial-  and  boundary- value  problem 
which  one  would  wish  to  solve,  would  be  well  posed.   In  particu- 
lar, if  in  any  sense  the  temperature  field  determines  the  wind 
field  and  if  a  state  of  quasi-geostrophic  equilibrium,  meaning  a 
state  of  minimum  amplitude  gravity-inertia  waves,  exists,  these 
two  facts  ought  to  reflect  an  instantaneous  property  of  the 
primitive  equations,  rather  than  an  integral  property  involving 
a  finite  time  span.   This  point  will  be  further  discussed  in 
Section  4.   Such  an  instantaneous  property  would  have  to  be 
expressed  mathematically  as  one  or  more  equations  involving  the 
velocity  components  and  their  space  derivatives,  but  not  their 
time  derivatives,  and  some  variables  of  state  (such  as  tempera- 
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ture  and  pressure),  which  are  assiimed  to  be  known  together  with 
any  number  of  their  space  and  time  derivatives. 

The  presence  of  as  many  of  the  time  derivatives  of  state 
variables  as  necessary  would  generalize  the  usual  forms  of  the 
balance  equation,  which  are  completely  free  of  time  derivatives, 
and  these  derivatives  would  be  all  that  is  needed  to  express  the 
thermodynamic  history  of  the  system.   These  time  derivatives 
could  be  computed  from  sufficiently  many  discrete  measurements 
by  using  adequate  numerical  techniques,  which  probably  would  have 
to  include  certain  space- averaging  procedures. 

If  the  number  of  independent  equations  thus  obtained  from 
the  primitive  equations  by  differentiation  and  elimination,  to- 
gether with  proper  side  conditions,  suffices  to  determine  the 
wind  field  at  a  given  instant,  then  we  will  have  solved  the 
completeness-of-data  problem,  assuming  adequate  knowledge  of  the 
thermodynamic  history.   If  not,  the  equation  or  equations 
obtained  will  still  express  a  compatibility  condition,  equiva- 
lent in  the  general  baroclinic  case  to  the  balance  equation. 
This  condition  will  then  be  the  only  limitation  imposed  by  the 
equations  themselves  on  the  initial  data  and  thus  the  only 
"dynamic"  check  on  the  accuracy  of  these  data. 

It  is  hoped  that  the  following  analysis  will  clarify  this 
point  of  view  and  indicate  some  of  its  possibilities  and  limita- 
tions. 


2.  Compatibility  conditions  for  the  shallow-fluid  equations 

As  an  example  of  the  ideas  involved  in  the  approach  we 
want  to  discuss,  we  consider  the  shallow- fluid  equations  for  a 
rotating  Cartesian  x, y-coordinate  system 

(la)  u^  "^^^x  '^^\  +  ^^-  fv  =  0   , 

( lb  )  v^  +  uv^  +  vVy  +  (t)y  +  f  u  =  0  , 

(IC)  (j)^  +U(t)^  +V(t.y +  (t)(u^  +  Vy)  =  0  , 

(t>  =  gh  . 

Here  u,  v  are  the  velocity  components  in  the  x,y  directions 
respectively,  h  is  the  height  of  the  free  surface,  g  the 
acceleration  of  gravity  and  f  the  Coriolis  parameter.   For  the 
sake  of  simplicity  we  assume  f  to  be  a  constant. 

First  we  will  study  the  linearization  of  this  system 
around  a  state  of  rest.   With  an  obvious  change  of  notation  we 
then  have 

(2a)  u^  +  (t)^-  fv  =  0  , 

(  2b  )  v^  +  ^y  +  f  u  =  0  , 

(2c)  (f)^  +0(u^+ Vy)  -  0  , 

(2d)  $  =  const. 

Here  $  is  the  equilibrium  value  of  the  geopotential  gh  of  the 
free  surface,  and  (|)  stands  now  for  gh  -  $,  the  deviation  of  the 
geopotential  from  $. 
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Differentiating  (2c)  with  respect  to  t,     (2a)  and  (2b)  with 
respect  to  x  and  y,  we  obtain 

(3a)  ^tx  +  ^xx-^^x  =   0  ' 

(3c)  *tt+^(^xt+V^  =  °  • 

Substituting  u  ,  and  v  ,  from  (3a),  (3b)  into  (3c)  yields 

(4)  ^f(^x-^y)-^(*xx+^yy^+^tt  =  ^  ' 

Using  our  standing  assumptions  that  (J)  and  as  many  of  its 
derivatives  as  necessary  are  known,  we  obtain  from  (4)  and  (2c) 
the  following  system  of  two  first-order  partial  differential 
equations  for  u,  v 

^x  +  V  ^  "^t/"^  ' 
(5) 

Uy-v^  =  (cf)^^-OA(j))/'Df  , 

where  A   is   the    two-dimensional   Laplacian   operator, 

Ai  =   i      +  d) 
^        ^xx      ^yy 

Thus  we  see  that  in  this  very  simple  case  it  is  possible  to 
eliminate  the  time  derivatives  u,  ,  v,  and  obtain  the  instantane- 
ous compatibility  conditions  (5)» 

System  (5)  for  the  functions  u  and  v  is  elliptic,  i.e.,  it 
has  no  real  characteristic.  In  fact  it  is  just  a  set  of  inhomo- 
geneous  Cauchy-Riemann  equations  for  the  functions  v,  u,  and  by 
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cross-differentiation  these  equations  lead  to  a  Poisson  equation 
for  either  u  or  v.   A  well-posed  problem  for  (5)  would  be  the 
Dirichlet  problem.   This  would  mean  prescribing  u  say  on  a 
closed  contour  ^D  and  v  at  some  point  on  the  contour  or  in  its 
interior,  D. 

It  is  well  known  (Morel  et  al.,  1971^  Williamson  and 
Dickinson,  1972)  that  the  system  (2)  has  three  independent  plane- 
wave  solutions,  one  corresponding  to  slow  Rossby  waves,  the 
other  two  to  inertia-gravity  waves  propagating  in  opposite  direc- 
tions.  In  the  case  at  hand,  where  the  unperturbed  velocity  is 
zero,  the  Rossby  mode  is  stationary  and  the  inertia-gravity  waves 
have  phase-velocity 

c  =  ±  (k^O  +  f^)^/2/k  ,    k^  =  k^  +k^  , 

k  =  (kT,k„)  being  the  wave  vector.   Any  solution  of  (2)  can  be 
represented  by  a  series  expansion  in  these  plane  waves. 
Let  a  quantity  X  be  decomposed  as 


X  =  X  +  X 


where  (~)  is  a  time  average  and  hence  X,  being  stationary,  is 
representable  in  terms  of  the  Rossby  mode,  whereas  X'  stands  for 
the  inertia-gravity-mode  part  of  the  representation.   Then 

(2')  fu  =  -^y  ,    ^~  =\ 

is  a  geostrophically  balanced  stationary  solution,  satisfying 
eqs.  (2)  as  well  as  system  (5)  with  B/^t  =  0. 
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The  most  general  stationary  solution  of  (2)  will  satisfy 


(6a)  u^+Vy  =  0  , 


(6b)  V  - u  =  Ay 


where  )(;  =  ^/f.      The  solution  of  (6a)  is 


(7)  u  =  -i/^   ,        V  =  i/^   , 


with  Tp   an  arbitrary,  twice  continuously  differentiable  function. 
By   (6b)  i;^  has  to  satisfy  the  equation 

(8)  A(^-  X)  =  0  . 

Hence  if 

iJ  =  X  on  ^D  , 

then 

^  =  X     in  D  , 

where  D  is  some  domain  and  ^D  its  boundary.   Similarly,  if 

h  th   =   "b  y        on  ^D  , 
XV  n'^ 


where 


\   =  n.V 


and  n  is  the  unit  normal  to  ^D,  V  the  gradient,  then 
(*)  ^  =  X  +consto    in  D  . 

But  we  have  from  (7)  that 


^n^  =  n. (v, -u) 
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Therefore  the  boundary  condition 


(9)  n. (v, -u)  =  S^x   on  ^^   ' 


and  in  particular 

(10)  u  =  -x^  ,    V  =  Xy    on  SD  , 

implies,  by  {*),    that 

u  =  -x^  ,    V  =  Xy   in  D  . 

Hence  we  can  conclude  by  a  continuity  argument  that  the  only 
case  in  which  (9)  can  hold  is  that  (lO)  be  satisfied.   An 
obvious  special  instance  of  the  Neumann  problem  above  is  that  in 
which  the  boundary  SD  is  a  level  line  of  the  height  field, 

BD:   X  ~   const.,    h  =  const., 

in  which  case 

u  =  -Xy  ^    ^  =  X^        in  D  , 

provided  ^D  is  also  a  streamline  of  the  velocity  field  (u,v)  and 

u  +  V  =  ( S  X )    on  SD  . 

Generally  speaking  we  can  state  that  the  time-averaged 
solutions  (u,v,(j))  of  (2)  will  be  geostrophic,  i.e.,  satisfy  (2')^ 
in  a  domain  D  if  and  only  if  they  are  "geostrophic"  on  ^D. 
Moreover,  because  of  the  maximum  principle  for  the  Laplace  equa- 
tion, deviations  e  from  geostrophicity  in  the  domain  D, 
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e  =  max  {  |  f  u  +  ({)y  | ,  |  f  v  -  c})^  | }  , 

will  be  bounded  by  the  values  of  e"  on  SD. 

Thus  we  see  that  in  the  simple  flow  described  by  (2)  we 
can  solve  the  completeness  problem  in  initialization  by  deter- 
mining u,  V  from  the  compatibility  conditions  or  "generalized 
balance  equations"  (5)?  given  ^,    (j),  and  <{'+■)- • 

Moreover^  a  nearly  geostrophic  wind  field  can  be  obtained 
from  a  time-dependent  geopotential  (j)  in  the  following  manner. 
First  expand  the  measured  ^   into  plane  waves,  as  indicated 
before.   Then  suppress  the  inertia-gravity  waves,  i.e.,  set  their 
coefficients  to  zero,  and  call  the  new  function  thus  obtained  (f)'. 
Clearly  (()!  =  0,  (j),  ,  =  0,  so  that  instead  of  (5)  we  will  have  for 
u,  V  a  system  of  the  form  (6),  with  )(  =  <{)'/£"•   0^^  "the  boundary, 
however,  the  departure  e  from  geostrophicity, 

e  =  max  (  |  f  u  +  (jj^  | ,  |  f  v  -  4^  | }  , 

will  in  general  not  be  zero.   But  we  have  shown  that  for  system 
(6)  the  error  e  in  the  interior  of  D  will  be  bounded  by  the  value 
of  e  on  the  boundary  ^D.   Hence  it  suffices  to  adequately  modify 
u,  V  on  ^D  and/or  (j)'  on  SD  (and  near  ^D,  to  preserve  continuity) 
in  order  to  keep  e  throughout  D  within  prescribed  limits.   This 
procedure  is  in  a  sense  comparable  to  that  used  in  the  analysis 
of  Williamson  and  Dickinson  (1972). 

We  are  now  in  a  position  to  return  to  the  more  general 
case  (1)  and  ask  whether  the  complete  analysis  performed  for  (2) 
generalizes.   Unfortunately  the  answer  seems  to  be  that  it  does 
not. 
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To  obtain  compatibility  conditions  similar  to  (5)  for  the 
system  (l)  we  would  have  again  to  eliminate  u,,  v,  from  it  by 
differentiating  (ic)  with  respect  to  t.   Before  doing  this  let 
us  simplify  (Ic)  by  introducing 

cp  =  log  (j)  . 

Then  cp  and  its  derivatives  will  be  known,  since  ^   and  its  deriva- 
tives are  known  by  our  standing  assumption.   Dividing  now  ( Ic ) 
by  ({)  >  0,  it  becomes 

(11)  u  +v  +  cp  u  +  cp  v+cp,  =  0  , 

^       '  X   y  ^x   ^y   ^t     ' 

where  cp  ,  cp  ,  cp,  are  known.   Differentiating  (11)  with  respect 
X   y    t 

to  t,  (la)  and  (lb)  with  respect  to  x  and  y,  we  obtain 

(12a)       u,   +UU   +VU   +u+uv-fv+(l)   =0. 
^         '  tx    XX    xy   X   y  X    X   ^xx 

( 12b )       v^  +  uv   +  vv   +  u  V  +  V  +  f  u  +  i    =  0  , 
ty    xy    yy   y  X   y    y   ^yy 

( 12c  )       u^y  +  Vy^  +  cp^u^  +  cpyV^  +  cp^^u  +  cpy^v  +  cp^^  =  0  . 

Now  we  can  substitute  u, ,  v,  into  (I2c)  from  (la),  (lb) 
and  u  , ,  V  ,  from  (I2a),  (I2b)  to  get  a  new  equation  for  u,  v 
containing  no  time  derivatives  of  u,  v: 

2  2 

( 13  )    uu   +  vu   +  uv   +  vv   +  u  +  2u  V  +  V  +  f  ( u  -  V  ) 
^  ^     XX    xy    xy    yy   x    y  x   y   '  y   x ' 

+  cp  ( uu  +  vu  +6     -  f V )  +  cp  ( uv  +  vv  +6     +  f u ) 
^x^   X    y   ^x     ^   ^y^   X    y   y 

-  ^xt^  -  ^yt^  +  ^xx  +  *yy  "  '^'tt  =  °  * 
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Together  with  (11)  this  would  then  yield  a  set  of  compatibility 
conditions  similar  to  (5)o   We  have  however  to  verify  whether 
(11)  and  (13)  are  indeed  independent  equations.   Since  (13)  is  a 
second-order  equation,  the  standard  procedure  to  test  for  inde- 
pendence is  to  introduce  new  dependent  variables  in  order  to 
convert  (11)  and  (13)  into  a  first-order  system. 
Letting 

Uy  =  q  ,    v^  =  r  , 

one  such  equivalent  system  for  the  six  unknowns  u,  v,  p,  q,  r,    s 
would  be 

u^  =  P  , 

Vy  =  S  , 


(14) 


p  -  q  =  0  , 

s  -  r  =  0  , 
X   y 

2         2 
up  +  vp  +  us  +VS  +p  +  2qr  +  r  +f(q-r) 
■^x    y    X    y  \  -1   / 

+  cp^(up  +  vq  -  fv)  +cp  (ur  +  vs  +  fu)  -  cp^^u  -  9y^v 

=  CD ,  ,  -d)   -i   -(pi-mi  ) 
^tt  ^xx  ^yy  ^x^x   y^y 

where  we  differentiated  (11)  with  respect  to  x  to  obtain  the 
third  equation. 

A  curve  is  called  characteristic  for  the  system  (l4)  if 
a  linear  combination  of  the  equations  of  ( l4 )  contains  only 
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interior  differential  operators  for  that  curve,  i.e.,  directional 
derivatives  along  that  curve.   Define 


to  be  the  components  of  the  normal  to  a  characteristic  curve 

^(x,y)  =  const. 


of  the  system  (l4),  and  let 


T  =  u|  +  vri 


Then  the  characteristics  of  the  system  (l4)  will  be  given  by 


/; 


det 


1 


-i 


\ 


X 


=   0 


) 


(e.g.,  Courant  and  Hilbert,  1962).   Developing  this  determinant 
yields 


(15) 


I  Tl(^TlT  -  ^Tjt)  =     0 


Therefore,  every  curve  is  a  characteristic  curve.   Hence  the 
matrix  differential  operator  on  the  vector  (u, v, p, q, r, s )  defined 
by  (l4)  is  an  interior  operator  on  every  possible  curve  of  the 
x,y-plane,  in  other  words  the  equations  in  system  (l4),  and 
therefore  equations  (11)  and  (I3)  cannot  be  independent. 
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Comparing  the  manipulations  which  led  to  equation  (13)  to 
those  which  led  to  equation  {^),    we  could  have  anticipated  this 
result.   Indeed,  in  (12c)  we  had  to  use  both  (la),  (lb)  in  order 
to  eliminate  u,,  v,  and  (12a),  (I2b)  in  order  to  eliminate  u  ,, 
V  , .   But  (12a),  (12b)  were  derived  from  (la),  (lb)  by  differen- 
tiation  and  thus  we  can  say  that  we  had  to  use  (la),  (lb)  twice 
in  order  to  obtain  (13).   On  the  other  hand  only  one  use  of  (2a), 
(2b)  was  necessary  to  obtain  (4).   The  obvious  analogy  with  a 
similar  situation  in  handling  linear  algebraic  systems  is  the 
intuitive  content  of  our  result  (15)« 

To  understand  the  deeper  reason  for  the  failure  of  the 
procedure  which  was  successful  in  the  case  of  system  (2)  we  will 
have  to  consider  the  general  case  of  the  primitive  equations. 
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3 .  Compatibility  conditions  for  the  primitive  equations 

For  the  sake  of  simplicity  we  shall  discuss  the  primitive 
equations  in  x,y, p-coordinates,  but  the  same  results  can  easily 
be  shown  to  hold  in  x,y, z-  and  x,y, 0-coordinates,  as  they 
obviously  should.   We  take  the  equations  in  the  form 

(l6a)  u^+Vy+oip  =   0   , 

( 16b  )  U ,  +  UU   +  VU   +  ODU   -  f  V  =  -  (i)    , 

^         '  t    X    y    p        ^x  ' 

( 16c  )  v^  +  uv^  +  vVy  +  a^Vp  +  f  u  =  -  (t)y  , 

(l6d)  e.+ue+YG     +0^6  =      0   , 

u     X     y     p 

where  c()  is  the  geopotential  of  an  isobaric  surface,  6  is  poten- 
tial temperature  and  00  =  dp/dt  is  the  vertical  velocity  in  this 
coordinate  system. 

Note  that  whereas  this  system  is  totally  hyperbolic,  i.e., 
it  has  four  families  of  characteristics  (three  of  which  coincide, 
V.  discussion  of  characteristics  and  their  determination  infra) 
the  three-dimensional  hyperplane  t=0  in  the  four-space  (x,y,p,t) 
is  not  free.    The  plane  t  ==  0  is  however  free  for  the  primitive 
equations  in  x,y, z-  and  x,y, a-coordinates  (where  a  time  deriva- 
tive appears  also  in  the  continuity  equation).   Hence  a  Goursat 
problem  (John,  1971,  p.  191)  with  part  of  the  data  prescribed  on 


The  Cauchy  problem  for  a  first-order  hyperbolic  system  consists 
in  finding  a  solution  of  the  system  with  values  of  the  unknown 
variables  prescribed  on  some  hypersurface,   A  surface  is  then 
called  free  if  the  Cauchy  problem  with  data  on  it  has  a  solution 
"in  the  small"  (i.e.,  at  least  in  some  neighborhood  of  the 
surface),  characteristic  otherwise  (John,  1971,  P«  5^  ff.).   In 
particular  a  Cauchy  problem  with  data  on  t  =  0  is  called  a  pure 
initial-value  problem. 
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t  =  0,  part  on  p  =  0  is  adequate  for  our  system,  rather  than  a 
pure  initial-value  problem. 

Again  we  assume  that  (f)  is  known,  together  with  all  its 
necessary  derivatives,  from  the  hydrostatic  equation 

(l)p  =  -RT/p  , 

where  R  is  the  gas  constant.   We  note  in  passing  that  this 
discussion  will  be  valid  no   matter  what  known  right-hand  sides 
equations  (l6)  might  have. 

We  want  to  find  a  set  of  equations  in  u,  v,  ca,    (j),  9   con- 
taining no  time  derivatives  of  u,  v,  cd,  or,  as  we  said  before, 
a  set  of  compatibility  conditions  the  wind  field  has  to  satisfy 
instantaneously.   Since  no  equation  for  ao^  exists  in  our  system, 
it  is  natural  to  use  (l6d)  to  eliminate  oj, 

(17)  ^  =  -  -^  (^x^'^  V^^t^  • 

P        ^ 

Differentiating  this  with  respect  to  p  and  substituting  into 
(l6a)  we  obtain 

(18)  u  +0^^^u  +v  +e^^K    +e^^K  +  e''^K  =  -6^^^  , 

^   ^      X       p   y       p 
where 

e^i)  =  -e  /e    ,  e^^^   =  _e  /e  , 

x'  p  '  y^  p 

0^^^  =  e  e  /e^  -  e    /e    ,    e^^^  =  e  u    /e^  -  e  /e    , 
X  pp^  p   xp^  p  '         y  pp^  p   y^  p 

t  pp^  p   tp  p 

are  known  according  to  our  standing  assumption,  since  T  and  p  are, 
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To  eliminate  the  time  derivatives  of  u,  v  between  (18)  and 
(l6b),  (l6c),  we  differentiate  (18)  with  respect  to  t,  and  (l6b), 
(l6c)  with  respect  to  p  and  with  respect  to  x,  y,  to  yield 

(18.)       u^,+e(i^Up,+Vy^  +  e(2)^^^^,(3)^^^e(4)^^ 

+  0(1)^   +e|2)^   +ep)^  +  ei^)v  =  -e|5)  ^ 
tptpt  t  t-* 

(I9a)        u,+uu      +VU      +a)u       +u+uv+ucD-fv     =  -^        , 
^         '  xt  XX  xy  xp        X       y  X        p  X  x  ^xx 

(  IQt)  )  U     ,    +  UU         +  VU  +  (JOU         +UU      +UV      +UCD      -fV        =-(j)  , 

^         ^  pt  xp  yp  PP        X  p        y  p        p   p  p  ^xp 

( IQc )        v,+uv       +VV       +a)V       +uv+v+va)+fu     =  -^        , 
^  yt  xy  yy  yp        y  x       y        p  y  y  ^yy 

( IQd )        v^+uv       +VV       +CUV       +uv+vv+vco+fu     =  -h 
^         ^  pt  xp  yp  PP        p   X       y  p        p  p  p  ^yp 

Substituting  u,  u,  ^   f    ^t  ^^'^'^    (19)  and  u,,  v  from  (l6b), 
(l6c)  into  (18')  we  obtain 

(20)    uu   +  VU   +a)U   +  uv   +  vv   +  odv 

XX    xy    xp    xy    yy    yp 

+  e'"^Muu   +vu   +(inu   )+e^^^(uv   +  vv  ^+cov^^) 
xp    yp    pp^       ^   xp    yp    pp' 

2  2 

+   u    +UV     +ua)    +UV     +v     +vtD 
xyxpxyxypy 

+   e^Muu     +UV     +uaj)+e^Muv     +VV     +V(d) 
^xp        yp        pp^  ^px       yp        pp' 

+   e^-^Muu    +vu    +(iDU    )+0^    Muv    +  vv    +(m) 
^      X  y  p^  ^      X  y  p' 

+  f(u   -V  )+f(0^^^u   -e^^^v  )-e|^^u   -e['^K 

^    Y        x'         ^  p  p^         t        p        t        p 


fie^^K-e^^K).ei^K-ei^^ 


V 


_^     .^     _0(i)^     .0(2)^     _0(3)^    _e(^)^   +e(5)   . 

^xx      ^yy  ^xp  ^yp  ^x  ^y        t 
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Here  again  cd,  cd   ,    (X)  ,  ca     can  be  eliminated  with  the,  aid  of  (17), 
^   y   P 

and  only  first-order  space  derivatives  of  u^,  v  appear  in  the 
process.   Thus  we  can  view  (20)  as  a  second-order  equation  in 
u,  V. 

To  test  the  independence  of  (18)  and  (20)  it  is  convenient 
to  introduce 

u  =  m  ,    u  =  r  , 

v  =  n  ,    V  =  t  . 
X     '     p     ^ 

and  differentiate  (18)  with  respect  to  p.   Then  we  obtain  the 
following  first-order  system 


u  =  Q  3 


(21) 


V   =  s  , 

y 


^x  -  S^  =  0  ^ 


n  -  s   =  0  , 
y   X 


q  -  r  =  0  , 

^p    X      ' 


s  -  t   =  0  , 

p   y 


uq 


(D...  .fl(i) 


+  vq  +  (oo  +  0^^u)q^  +  e^  ^(vr  +  cor  ) 
X    ^   '         ^^P       ^   y    P^ 

+  US  +VS  +(aj  +  e^^^v)s  +e^^^(ut  +cot  )  +  .o.=0  , 
X    y   ^         '  p       ^   X    p  ^         ' 

r    +  e^-^K    +  t  +e^^K    +  ,  ..  =0  , 
X       p   y       p         ' 


where  the  dots  stand  for  terms  which  contain  only  the  unknowns 


23 


U  =  (U;,  v,m,  n,  q,  r,  s,  t )  in  non-differentiated  form  or  not  at  all. 
Thus  system  (21)  can  be  written  as 

A,U  +A^U  +A^U  +B  =  0  , 
1  X   2  y  ^  p 

where  A,,  k^,    A^  are  matrices  depending  on  U^  but  not  on  its 
derivatives  and  B  is  a  column  vector  with  the  same  property. 
Let  a  surface  in  x,y,p- space  be  given  by  an  equation 

ip{y^,Y  }V)   =   const, 
and  denote  by 

^  =.  ^^  ,   n  =  ^y  ,   ^  =  ^p 

the   components   of   its   normal.      Consider  the   matrix 
A   =  A-^^  +A2Ti+A^C    = 

C  -I 

T+(a)+e^^^u)C      e^^^VTi+ojC)      T+(a3+e^^V)C      0^^''(u^+a)C 


where  we  used  a  slight  reordering  of  the  equations  of  (21)  (the 
equation  which  was  considered  the  fifth  one  in  the  first  writing 
being  now  last)  and  where  t  =  u|  +vrio 
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By  the  theory  of  characteristics  the  quasi-linear  system 

(21)  is  elliptic  at  a  point  (x,y, p)  for  a  certain  solution 
U  =  U(x, y,p)  if  the  equation 

(22)  det  A(x,y,p;U;4,ri,C)  =  0 

has  no  real  solution  (^,ri,C)  there  and  is  totally  hyperbolic  if 
eight  distinct  real  solutions  {^,T],t^).,    i  =  1,2,... ,8,  exist. 
Any  such  "solution"  of  (22)  is  a  first-order  partial  differential 
equation  for  a  function  il/ .      The  solution  of  this  PDE  yields  a 
family  of  characteristic  surfaces 

^(x^y^p)  -   const. 

of  the  system  (21).   The  equations  of  (21)  cannot  be  independent 
of  each  other  if  (22)  is  identically  satisfied,  irrespective  of 
the  values  of  the  variables  on  which  A  depends,  as  was  already 
explained  for  the  shallow-fluid  equations.   Indeed,  for  (22)  to 
be  identically  satisfied,  it  is  necessary  that  the  matrix 
differential  operator 

C  =  A-^V^x  +  A^V^y  +  A^V^P 

be  an  interior  operator  at  any  point  of  any  surface  in  space  for 
any  solution  U  of  (21).  In  other  words,  a  solution  U  of  some  of 
the  equations  in  system  (21)  must  satisfy  the  whole  system 

C(U)U  +  B(U)  =  0  . 

Developing  the  determinant  in  (22)  by  minors  yields 


25 


-  e^2^[u|+VTi  +  (a3+0^^^u)C]^n  =  0  . 

Thus  the  equations  of  system  (21)  are  not  independent.   Hence 
the  same  holds  for  equations  (18)  and  (20),  to  which  system  (21) 
is  equivalent,  given  suitable  side  conditions. 

We  remark  moreover  that  neither  is  a  full  set  of  compati- 
bility conditions  possible  for  the  complete  Euler  equations,  in 
which  the  equation  for  the  vertical  velocity  cannot  be  reduced 
to  the  hydrostatic  equation.   Indeed,  let  us  assume  that  somehow 
all  variables  of  state,  i.eo,  T,  p,  0   are  known  together  with 
the  necessary  derivatives.   Even  then,  differentiating  the  energy 
equation  equivalent  in  x, y, z-coordinates  to  (l6d)  with  respect  to 
t  and  substituting  u,  ,  v,  ,  w,  from  the  three  momentum  equations 
into  the  equation  so  obtained  will  lead  to  one  compatibility 
condition.   The  continuity  equation  will  then  give  just  one  addi- 
tional condition,  and  it  is  verified  by  inspection  that  no  other 
independent  equation  in  which  time  derivatives  of  (u, v,w)  do  not 
appear  can  be  derived.   Thus  we  have  again  an  underdetermined  set 
of  two  equations  for  the  instantaneous  determination  of  the  three 
unknowns  u,  v,  w. 
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In  this  context  it  is  worthwhile  to  note  that  in  Mesinger''s 
(1972)  iterative  updating  experiments  the  use  of  the  first  time 
derivative  of  the  state  variable  brought  a  significant  improve- 
ment in  the  adjustment  of  the  wind  field.   The  use  of  the  second 
time  derivative  however  did  not.   This  reflects  the  fact,  which 
should  be  clear  in  the  light  of  our  discussion,  that  no  addi- 
tional independent  compatibility  conditions  can  be  obtained  by 
further  differentiation  with  respect  to  time  of  the  energy  or 
continuity  equations. 

Obviously,  to  consider  the  Navier-Stokes  equations,  or 
other  equations  with  forcing  and  dissipation  depending  on  the 
wind  field,  would  make  matters  only  worse.   This  seems  to  exhaust 
the  possibilities. 

We  are  therefore  forced  to  conclude  that  there  is  no 
adequate  generalization  of  exact  geostrophic  balance  to  atmo- 
spheric phenomena  governed  by  the  primitive  Euler  or  Navier- 
Stokes  equations,  at  least  not  in  the  sense  that  the  wind  field 
is  entirely  and  instantaneously  determined  by  the  thermodynamic 
state  of  the  atmosphere  and  its  past  history  as  reflected  in  the 
time  derivatives  of  the  state  variables. 
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4.  Discussion  and  conclusions 

The  analysis  above  should  make  it  clear  why  the  mathemati- 
cal argument  mentioned  in  the  Introduction  and  used  by  some 
authors  to  justify  the  updating  approach  to  the  initialization 
problem  may  not  be  appropriate.   Not  only  are  the  systems  of 
PDEs  under  consideration  non-linear,  which  in  general  precludes 
the  elimination  of  all  but  one  variable;  the  even  more  pertinent 
difficulty  is  that  such  elimination  does  not  lead  to  a  well- 
posed  problem  to  determine  the  wind  field  from  temperature  data. 

On  the  other  hand,  our  analysis  proves  rather  conclusively 
that  for  an  atmosphere  governed  by  the  primitive  equations  even 
a  generalized  geostrophic  balance,  in  the  sense  of  compatibility 
conditions  as  previously  defined,  can  be  only  an  approximate, 
rather  than  an  exact  mathematical  property. 

¥e  think  that  this  discussion  at  least  partly  explains  the 
difficulties  encountered  by  the  updating  and  four-dimensional 
data  assimilation  approaches  to  the  initialization  problem. 
Specifically  these  difficulties  are:   (i)  that  the  erros  in  the 
initial  data  cannot  be  reduced  to  zero  or  to  negligible  values, 
(ii)  that  in  forward  updating  the  asymptotic  value  of  the  error 
is  reached  after  times  comparable  to  the  predictability  time,  and 
(iii)  that  use  of  data  extraneous  to  the  numerical  model  used 
seems  to  lead  to  catastrophic  results.   The  partial  success  of 
these  approaches  has  probably  to  do  with  a  rather  complicated 
mechanism,  involving  the  forcing  and  dissipation  in  the  system 
of  equations  used,  which  we  cannot  discuss  here  and  which  is 
currently  under  investigation. 
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Our  analysis  on  the  other  hand  has  shown  that  instanta- 
neous compatibility  conditions  for  the  velocity  field  can  he 
derived,  but  that  their  number  (one  for  the  non- linear  shallow- 
fluid  equations,  two  for  the  primitive  equations)  is  in  general 
one  less  than  the  number  required  to  determine  the  velocity 
field  from  the  thermodynamic  fields  and  their  history.   This 
seems  to  justify  the  contention  that  it  is  necessary,  in  order 
to  solve  the  initialization  problem,  to  be  able  to  measure  at 
least  one  scalar  parameter  determining  the  wind  field.   This 
scalar  could  be  either  the  (horizontal)  wind  direction  or  the 
wind  speed. 

The  wind  direction  could  possibly  be  determined  by  cloud 
tracking  techniques  which  are  not  at  present  sufficiently  accu- 
rate in  order  to  determine  also  the  speed.   The  speed  on  the 
other  hand  could  perhaps  be  determined  with  sufficient  accuracy 
(concerning  the  present  GARP  requirements  see  for  instance 
Kasahara,  1972)  by  other  methods  which  are  still  easier  to 
implement  than  those  requiring  a  full  and  accurate  measurement 
of  two  or  three  velocity  components  separately. 

Our  analysis  may  be  extended  to  study  four-dimensional 
data  assimilation  in  the  following  way.   Assume  that  one  would 
wish  to  solve  a  Cauchy  problem  for  a  hyperbolic  system  of  equa- 
tions of  large-scale  dynamic  meteorology  (the  shallow-fluid  equa- 
tions or  the  primitive  equations)  different  from  the  pure 
initial-value  problem.   This  would  correspond  to  data  being 
available  on  some  arbitrary  hypersurface  S  in  (x,y, z,t)  four- 
space,  rather  than  on  the  hyperplane  t  =  0,  e.g.,  on  a  surface 
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spanned  by  the  trajectories  of  satellites.   If  the  data  were 
complete  and  the  given  surface  S  were  free  for  the  system  of 
equations  used,  as  well  as  being  a  space-like  surface  (John, 
1971,  p.  115  ff.,  Courant  and  Hilbert,  1962,  p.  589  ff . ),  then 
the  Cauchy  problem  with  respect  to  S  would  be  well  posed. 

If  the  data,  however,  are  not  complete,  as  is  to  be 
expected,  although  S  is  free  and  space-like,  we  have  again  with 
respect  to  S  an  "initialization  problem"  or  "incomplete  Cauchy 
problem".   This  can  now  be  handled  either  by  "updating  tech- 
niques", as  is  currently  done  in  various  numerical  experiments, 
or  else  by  the  analytical  methods  proposed  above.   Our  analysis 
of  the  initial-value  problem  seems  to  suggest  that  a  theoretical 
investigation  of  four-dimensional  data  assimilation  in  the  sense 
outlined  here  might  help  direct  the  n-umerical  experiments  in 
such  respects  as  optimal  number,  nature  and  location  of  data,  or 
assessment  of  the  completeness  and  relative  value  of  data  avail- 
able under  present  and  future  data  gathering  programs. 

This  subject,  as  well  as  the  numerical  implementation  and 
practical  evaluation  of  the  theoretical  methods  proposed  here 
for  the  initial- value  problem,  are  intended  to  be  covered  in 
future  work  by  the  author. 
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